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A low-dimensional model (LDM) for turbulent Rayleigh-Benard convection in a Cartesian cell 
with square domain, based on the Galerkin projection of the Boussinesq equations onto a finite 
set of empirical eigenfunctions, is presented. The empirical eigenfunctions are obtained from a 
joint Proper Orthogonal Decomposition (POD) of the velocity and temperature fields using the 
Snapshot Method on the basis of a direct numerical simulation (DNS). The resulting LDM is a 
quadratic inhomogeneous system of coupled ordinary differential equations which we use to describe 
the long-time temporal evolution of the large-scale mode amplitudes for a Rayleigh number of 10 5 
and a Prandtl number of 0.7. The truncation to a finite number of degrees of freedom, that does 
not exceed a number of 310 for the present, requires the additional implementation of an eddy 
viscosity-diffusivity to capture the missing dissipation of the small-scale modes. The magnitude 
. of this additional dissipation mechanism is determined by requiring statistical stationarity and a 

total dissipation that corresponds with the original DNS data. We compare the performance of 
two models, a constant so-called Heisenberg viscosity-diffusivity and a mode-dependent or modal 
one. The latter viscosity-diffusivity model turns out to reproduce the large-scale properties of the 
turbulent convection qualitatively well, even for a model with only a few hundred POD modes. 
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For most turbulent flows in nature and technology, it is impossible to resolve all relevant degrees of freedom. 
Systematic methods to derive models with a reduced number of degrees of freedom from the full set of nonlinear fluid 
£L| equations are thus necessary. Low-dimensional modeling of transient and turbulent flows using Galerkin projection 
onto the empirical basis functions which are obtained from a proper orthogonal decomposition (POD) is one such 
well established method^— POD and the development of Galerkin models based on POD modes has been applied 
■ to a number of fundamental hydrodynamic flow problems, including simple wall-bounded shear flows* 4 -"— flows over 
t*"*- \ cavities^ or in the wake of a cylinder fi£~— The development of low-dimensional models (LDM) based on POD modes 
i has also been extended in several directions such as to the balanced POD method- 4 !- 5 or to unsteady flow problems 15 
for which fast and slow flow modes are separated. Most of these cases have been studied for laminar or transitional 
, flows at lower or moderate Reynolds numbers. 

With increasing Reynolds number the flows become turbulent, the number of degrees of freedom grows rapidly 
and their nonlinear couplings are increasingly relevant. The truncation of the set of nonlinear ordinary differential 
equations (ODE) which follows from Galerkin projection introduces always a cut-off of these mode interactions and 
removes couplings between the degrees of freedom which are necessary for the transfer of kinetic energy from large to 
small scales. An additional dissipation mechanism has to be implemented in the low-dimensional model to account 
k> ' for the dominant dissipation by the truncated degrees of freedom. The particular way of truncation can then alter 
. the dynamics in the LDM. 

Several approaches to this problem have been suggested in the past. Aubry et alA used directly the energy transfer 
between resolved modes and unresolved modes at smaller scales to formulate a spectral closure in the truncated system 
of ODEs. Moehlis et aZ.— presented streamwise-invariant truncations for a plane Couette flow, and showed that very 
low-dimensional models with up to ten degrees of freedom can reproduce transient flow phenomena in low-Reynolds- 
number shear flows. They also found that the detailed behavior of their LDM depends in a subtle manner on the 
modes included and that a proper account of the symmetries of the system is crucial. Later Smith et ai£^ included 
streamwise variations in their model. They also introduced a linear damping term, but only when the particular 
POD mode expansion coefficient is significantly anti-correlated with its time derivative, da^ jdt. Cazemier et 
ai& constructed a LDM for driven cavity flows, consisting of the 80 most energetic POD modes computed from 700 
snapshots of a direct numerical simulation (DNS). To study the time evolution of the truncated ODE system, a 
slightly different linear damping term is introduced in their model. This term is calculated from the requirement that 
the energy of the ODE system is conserved in a statistically stationary sense. 

A few attempts to derive LDMs are reported for Rayleigh-Benard (RB) convection, despite being one of the most 
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comprehensively studied flows i 17 ' 18 Studies of RB convection in a finite box, based on the POD procedure, have 
been mostly done by Sirovich and co-workers i 1 ^— Sirovich and Park22i2i discussed the importance of the discrete 
symmetries describing the velocity-temperature fluctuations field. Deane and Sirovich 2 ^ made a parametric study of 
the POD mode spectra for small Rayleigh numbers Ra < 46000. 

Only recently, a snapshot method has been applied to turbulent RB convection in a closed cylindrical cell for 
Rayleigh numbers up to Ra = 10 s and cell aspect ratios between one half and threeJs In this work, emphasis was 
given to relating the first POD modes to the large-scale flow circulation which is always present in a closed turbulent 
convection cellJ^ The disentanglement of the temperature and velocity fields into POD modes allowed the authors to 
quantify the amount of heat which is transported by the particular POD modes through the convection cell. A change 
of the large-scale flow from a one-roll to a two- roll pattern, which is observed when the aspect ratio is increased beyond 
one at a fixed Rayleigh number, was in line with a decrease of transported heat by the primary mode compared to 
the secondary POD mode. 

As a correspondence of the few POD studies of RB convection, only a few works exist with an emphasis on developing 
a LDM by a Galerkin projection of the Boussinesq equations onto the most energetic POD modes. Tarman^ derived 
a model from POD modes which have been however separately extracted from the velocity and temperature fields. In 
a second work he proposed an algorithm which incorporates the lost dissipation due to truncation^ Besides the cutoff 
index based on the energy (mode index k < fc e ), a second index based on the dissipation (kd > k e ) was considered. 
The time dependence of the modes with indexes k e < k < kd was expressed as the quotient of the corresponding 
nonlinear and dissipation coefficient. No closed forms for the constant coefficients in the ODE system were however 
obtained in any of these works. 

In the present work, we want to extend these studies of RB convection in several directions. First, we construct a 
LDM for the evolution of the POD mode coefficients in the case of turbulent Rayleigh-Benard convection in a 
Cartesian cell with periodic side walls and isothermal free-slip top and bottom square planes. It is essential to use 
POD modes of the combined four-vector velocity-temperature field. 24 In this derivation, it turns out that a cubic term 
due to the interaction of the velocity with the mean temperature field (denoted as 9) becomes linear as a consequence 
of the orthogonality of the POD modes. The other terms which arise in the Galerkin projection are a linear production 
term p, a linear dissipation term e, a quadratic nonlinear term N, and a constant term em corresponding to the 
dissipation due to the mean temperature field. Second, our studies will extend previous work o 21 ' 22 i 26 in terms of the 
magnitude of the Rayleigh number of convection. A case with Ra ~ 10 5 is considered for which RB convection is 
turbulent and a DNS data record exists. Third, we are interested in the long-time behavior of the dynamics in the 
LDM. With a view to more complex convection flows in the future, we are seeking for the least set of POD modes 
that can reproduce characteristic dynamics of turbulent convection. 

A solution which includes the additional dissipation due to the neglected less energetic POD modes has to be 
considered by an additional eddy viscosity-diffusivity, r\ > 0. First, we present the so-called Heisenberg model with a 
constant r\ which exerts the same fraction of dissipation on all POD modes. As will be shown, this closure requires at 
least a minimum number of degrees of freedom, in particular with respect to the vertical direction, for a qualitatively 
correct description of the flow. As a consequence, two LDMs with 210 and 310 degrees of freedom, respectively, are 
chosen. They are taken from a set of 15708 modes (see Sec. Ill B). As will be seen, this model fails to reproduce 
the large-scale evolution of convection. For the larger of the two sets of modes, the model relaxes to a statistically 
stationary state which contains too much energy. Second, we refine this model and include a mode-dependent (or 
modal) eddy viscosity-diffusivity. The magnitude of both eddy viscosity-diffusivity contributions has to be estimated. 
In order to do so, we will follow a procedure that has been suggested by Cazemier et al&. The second model yields 
much more realistic large-scale variations of the most energetic modes, also reproducing with reasonable accuracy the 
energy spectrum and the turbulence statistics. Therefore, a significant part of the present work discusses the impact 
of both types of eddy viscosity-diffusivity rj on the dynamics of the LDM with different number of degrees of freedom 
and how it compares to DNS. 

The outline of the paper is as follows. The equations of motion, the basic idea of POD - in particular for the 
method of snapshots - is discussed in the next section. The construction of the LDM by Galerkin projection onto 
POD modes of RB convection follows in Sec. III. In this section, the results of the time integration of the LDM with 
both eddy viscosity-diffusivity schemes, and the agreement with the DNS are also discussed. We conclude with a 
summary and give an outlook. 
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II. METHODS 
A. Equations of motion and numerical scheme 

Turbulent Rayleigh-Benard convection is governed by the Boussinesq equations. They are brought into a dimen- 
sionless form by rescaling with the domain height l Zl the diffusive time scale = l z /n with n being the thermal 
diffusivity, the temperature difference AT — Qbottom — ©top > and follow to 

Vu = (f) 
+ (u • V)u = -Vp + PrV 2 u + RaPrTe z (2) 

at 

dT 

— + (u- V)T = w + V 2 T (3) 

where u(x, t) is the velocity field, T(x, t) the departure from the linear conduction temperature profile, and p(x, t) 
is the kinematic pressure. Dimensionless parameters are the Rayleigh number Ra — gaATl z 3 /(vn) and the Prandtl 
number Pr — v/k. Besides diffusivity k, they contain the kinematic viscosity v, the gravitational acceleration g, and 
the thermal expansion coefficient a. Note that the total temperature field is given in our notation by (see also Ref<2£) 

AT 

&(x,t) = Q bottom - — z + T(x,i). (4) 

The vector e z is the direction in which buoyancy and gravity work and in which the mean temperature gradient is 
established. The dimensions of the cell are l x = ly = 4-7T, l z = ir, where from now on x, y are the horizontal and z the 
vertical dimensionless coordinates. The aspect ratio is fixed to l x /l z = l y /l z = 4. For convenience, the origin of the 
coordinate system is in the center of the cell. Therefore, x £ [— L x /2, L x /2], y € [— L y /2, L y /2], and z E [— f/2,f/2] 
with L x = l x /l z and L y = l y /l z . The x-, y- and z-components of the velocity field will be denoted by u, v and w, 
respectively. The boundary conditions are periodic in x and y, and free-slip in z. This means that at the hot bottom 
plane at z/l z = —1/2 and the cold top plane at z/l z = 1/2 the following conditions hold: 

w = T= — = — = 0. (5) 

oz oz 

For the present boundary conditions, the flow can be decomposed in 

m(x, t) = (u(z)} + w'(x, t) (6a) 

v(x,t) = {v(z))+v'(x,t) (6b) 

w(x, t) = (w(z)) + w'(x, i) (6c) 

T(x,t) = (T(z)) + 9(x,t) (6d) 

where the mean components, e.g. (u(z)} are ensemble averages obtained by averaging over the horizontal x-y plane 
and time, i.e., a sequence of N' T statistically independent snapshots. The ensemble average for the velocity component 
u is thus given by 
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r W2 ,L y /2 \ 

( u i z )) — ttt I T~i~ I / u(x,y,z,t n )dydx . (7) 

T n=1 \ L ^ L y J -L v /2 J 

Figure [1] shows the mean vertical profiles of and T together with the linear thermal conduction profile. Considering 
the fact that for our problem, due to symmetry considerations, (u) = (v) = (w) = 0, we can take the four-vector field 

v(x,t) = (u,v,w, 9) (8) 

for the POD analysis and LDM derivation. Our analysis is in the statistically stationary regime of convective turbu- 
lence. The ensemble average of ([3]) yields then an expression linking the mean and the fluctuating components of the 
flow in the form 

d 2 (T) _ d^w9)_ 

dz 2 ~ dz 1 ' 
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FIG. 1: Mean vertical profiles (T(z)), (0(z)) and the linear conduction profile (see Eq. Q). Data are obtained from the DNS 
data record. 



Following Sirovich et a/.— temperatures and velocities are rescaled with 6 C = ^/2Nu 3 / Ra Pr and u c = \J Ra Pr/2Nu, 
respectively. Here Nu is the Nusselt number. The quantity C is obtained by demanding that the turbulent heat flux 
in the center of the cell must be equal to that due to diffusion at the boundary. 

The highly-resolved data record is obtained by a pseudo-spectral DNS which uses fast Fourier transformations . 27 ! 28 
Time stepping is done by a second-order Runge Kutta scheme. For most of the work, we consider a data set with 
N!p = 320 full three-dimensional turbulence snapshots which are separated by four convective time units t conv — 
\fl z / gaAT from each other. The computational grid consists of N x x N y x N z = 256 x 256 x 65 points in x-, y-, 
and ^-directions, respectively. The spectral resolution is given by k max r\K = 5.6. Here, k max = \/2N z /3 and t\k the 
Kolmogorov dissipation length. The Rayleigh number is Ra — 1.03 x 10 5 and the Prandtl number Pr — 0.7. 



B. Proper Orthogonal Decomposition (POD) 

The POD is a model reduction technique that extracts the most energetic modes from a set of realizations or 
snapshots of the flow. These POD modes are used as a basis for Galerkin projections of the full set of nonlinear 
equations thus reducing the infinite-dimensional space of solutions to a finite-dimensional system* 2 -^ The two-point 
correlation tensor or covariance matrix of the four- vector field is defined by 

X m „(x,x') = (v m (x,t)v* n (x',t)) t (10) 

where the asterisk denotes the complex conjugate, (-)t the time average and m, n = 1,2,3,4. For Rayleigh-Benard 
convection in Cartesian domains with two homogeneous (invariant with respect to translations) directions Eq. (p~0|) 
takes the forrr>2£ 

KT mn (x,x') = K mn (x - x',y - y',z) (11) 
For a kernel {jXTJ) , the eigenfunctions have the form 

m (,, y, z) = exp + *Jpy) (12) 

V / L X Ly \ L X Ly J 

where n Xl n y are integers for the a;— and y— directions, respectively. The superscript (p) denotes a particular POD 
mode. The determination of 4> follows then from 

4 „l/2 

V/ K mn (n x ,n y ;z,z')^) nxtn (z')dz' = (z) (13) 

where K mn is the Fourier transform of K mn with respect to the homogeneous directions x and y. The kernel K mn 
is calculated from the numerical data set by first taking the discrete Fourier transform of each realization in the 
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dxdy (14) 



horizontal plane, 

1 f Lj2 [ Lv/2 , N r f2mn x x 2nin v y\ 
F m (n x ,n y ;z,t) = — — / / v m (x, y, z,i)exp - — 1 j 

L x L y J_l x /2 J-L y /2 I \ L y J . 

and then averaging the correlation over the entire ensemble of data, 

K mn (n x ,n y ;z,z') = (F m (n x ,n y ; z,t)F*(n x ,n y ; z' ,t)) t . (15) 

Thus the kernel K mn is Hermitian, non-negative and on physical grounds square integrable, such that the existence of 
a complete set of vector eigenfunctions {4>m]n x ,n v (z)}p=i... given by (fl~3|) is assured. Complex conjugation of Eq. (p~3|) 



and use of Eq. (jl2j) implies that 

Due to the reality of the physical space fields n m (x, t) Eq. (TTJJ implies that 

F^(n x ,n y ;z,t) = F m (~n x , -n y ; z,t) . (17) 
The associated expansion of the velocity field w m (x, t) in terms of the modes is given as 

/ ,s V^V^ OnJ,n y (t) f 2irin x x 2mn v y\ { ) 
v m (x,t)=}^}^}^—==exp[— + —T- )^L,n v (z), (18) 

p n x n y \l L * L y \ ^ / 

and again reality of the four- vector field implies that cinl,n y {t) = a -«* -n W- The index (p) runs over the POD 
modes. The coefficients are calculated by the scalar product in L2(f2), 

ai nln y (t) = fol^). ( 19 ) 

4 ,1^/2 /.i v /2 ,1/2 

= V / / / {x,y,z)v m (x,y,z,t)dxdydz. 

m=1 J-L I /2J-L v /2J-l/2 

Next, the discrete Fourier transform of w m (x, t) is introduced together with the ansatz of <J> in (|12l) and the orthogo- 
nality of the complex exponentials. This gives 



m=l J ~ 1 



1/2 

CLnJ 2 )^^' n y ;z,t)dz. (20) 

/2 



C. The method of snapshots 

The snapshot method is one way to obtain the POD modes, particularly when the computational grid size becomes 
large. The time coordinate t in the equations above has to be substituted now by an index that runs over the sequence 
of snapshots. It is based on the fact that (TT51) is a degenerate kernel.— Consequently, an eigenfunction of the kernel 
K mn can be represented as 

0mL,™«( z ) = ^2^2^ ■ a {p) (n X7 n y ,l)j ■ F m (n X7 n y ;z,l) , 

1=1 7£S 
N T 

= ^a ( - p \n x ,n y ,k)F m (n x ,n y ;z,k) , (21) 
fe=l 

where 7 as an element of the symmetry group Q = {G, ZG} as discussed in AppendixlAl The explicit use of symmetries 
in the problem at hand enlarges the data record with originally N T snapshots to a total number of Nt — 16 
snapshots thus improving the convergence. For N' T = 320 DNS snapshots we thus end up with Nt = 5120 samples 
that can be used to evaluate the POD modes. Replacing the kernel in (fT3| and using Eq. (j2~Tj) results to 



L x L y ( f F*(n x ,n y ;z', k)Fj(n x , n y ; z' , m) dz'^] a {p) (n x , n v ,m) 

*T m=1 \J-l/2 J = 1 } 

\W a W(n x ,n v ,k), (22) 



6 



where k,m — 1,2,...,Nt represent any two snapshots (including all possible symmetries). Then (|22|) is the matrix 
problem which yields the eigenvalues A and eigenfunctions <j). It is clear that it determines just Nt of the empirical 
eigenfunctions for a fixed tupel (n x ,n y ). The eigenvalue Xn],n y of the Nt x Nt matrix is the total energy (kinetic 
energy plus temperature variance) of the pth POD mode for (n x , n y ). Recall that in the present case the four- velocity 
field is expanded into Fourier modes with respect to x and y which is characterized by wavenumbers n x and n y , 
respectively. This results theoretically in an infinite set of POD modes. 



III. RESULTS 

A. Galerkin Projection of the Boussinesq equations onto the POD modes 



Given the full set of nonlinear Boussinesq equations and the POD modes extracted by a snapshot method from the 
DNS data, we can proceed to derive the LDM. This requires first a Galerkin projection step. Using the dimensionless 
units introduced in section III A| Eqns. ([2]) and ([3]) can be rewritten together in a four- vector notation with respect to 
v, 



3=1 



< / > T\ Yl 7TT + 2Sl3 Nu2 V4+ 
+ 8ii(Pr - 1) ^ j 



0i4 ~T v 3 



A d 2 {T) 

2^ - 



i dx 2 



2^ Vj dx, 

j= i j 



(23) 



(24) 



where i = 1,2, 3, 4. Here, j = 1, 2, 3 correspond to the three spatial coordinates and the term 

s t = 25 a Nu 2 (T)-i|P , 
u c dxt 

is a source which drops out in the Galerkin projection procedure. This is due to the divergence-free nature of the 
POD basis functions and the fact that modes $g P Q (x, y, z) = 0, Vp, respectively. Now one takes the inner product of 
((221) with modes ^m} Px , Pv {x,y, z) and inserts the expansion v m (x,t) = ^2 n Y^n x a £«, (t)®m}n*,n v (x, y, z) from 
equation (jTSJ). Due to the orthogonality of the POD modes, the following infinite-dimensional ODE system follows 



■Jp) = 

Px ,Py 



n=l 



n=l 



N(p x ,Py,p) + ^2 ^(Px'Py> n ,P) a pl]p y + e {T)(Px,Py,p) 



(25) 



The terms on the right hand side of (I25|) correspond to the production (p), dissipation (e), nonlinear transfer (N), 
interaction with the mean flow (3), and the dissipation due to the mean temperature field (c(t>)j respectively. Closed 
forms for e and N have been obtained in RefsJ^ and^, but only for the velocity field. The general equations for the 
terms on the right hand side of (|25[) are as follows. The production term is given by 



p{Px,P y ,n,p) = 
and the dissipation term by 

e(Px,Py,n,p) 



1/2 



1/2 



2Nii 2 6 < - n) 6 (p> * + — 6 (n> 6 {p> * \dz 

ZlyU ^Mp^Py^Px^y + n P^PvPy^^P^Py 



k(p) 



Pr± 

i=i 



2irp x 



2-KPy 

Li, 



2irpy 









2-i 


,1/2 




/ 

J-l/2 




"1/2 , 


J 


-1/2 




2- 


,1/2 




/ 

J-l/2 




,1/2 
/ 


J 


-1/2 



dp)* 



(26) 



T,Px,Py ^r,p* ,Py 



dz+ 



(™)' (jS-P)*' 



dz 



dn) dp)* t , 

^■,P*,Py ( l > i\Px,Py aZ ^ 



in)' ,(p)*' , 
4;Px,P H ^i; P:c ,py U " : 



(27) 



7 



The nonlinear mode coupling term N is given by 

oo oo oo oo 



n—1 m— 1 m x — — oo m y — — oo 

where the coefficients are given by 



2mm y (n) (m) An) Am)' 



,1/2 


'2 


-1/2 









dn) Am) 

l;p x —m x ,p y —m y r'j-.m^ ,m y 



S, P /*' (29) 



The last two terms in Eq. (f2"5)l . 9 and £(t)> are calculated for a quasi-steady flow^. The term 9 is related to the 
ensemble average of temperature, u • V(T) = w^P-. However, from the ensemble average of ([3]) we obtain 

ffi(z) = u c {w9{z)) + ^(z= -1/2) (30) 
dz dz 

where -j^-{z = —1/2) = —u c j_^, 2 (wd) dz. In terms of the eigenvectors, 

oo oo oo oo 

<^> = tVEE E E (^w-V^c*))^.*^-*.,-*. ( 31 ) 



fc— 1 /— 1 fc^— — OO /Ci, — — OO 



since (u>#) is only a function of z. However (ai \ _t (*)) = ( a l^z, (*) a l fc W) = ^fcAt & an< ^ therefore 

1 OO OO OO oo 

- j^EE E E MS^U^-*..-*. ( 32 ) 

33 * fe=l i=l fc^— — OO fcy— — OO 

^ oo oo oo 

~T" E E E ^,fc y ^3;fc x ,fe H ^4;fe a; ,fc H ■ ( 33 ) 



L x L y - 

w k— 1 — — oo k y — — oo 

Consequently, the equation for the term 9 is 

9(p B ,ft„n,p) - -«= / 1/2 - dz) ^Lv^Lvy dz ( 34 ) 

J-l/2\ J -1/2 / 

where (u>0) is given by ([33]). For the low-dimensional description of free-shear-flow, Rajaee et al£^ keeps the original 
time dependence in (|33[) . just replacing Aj, fc = (|a^ fe | 2 ) and considering this as a running-time- average factor, so 

the resulting J2^Li ^(k x ,k y , n, k)a^ k term becomes cubic in the ODE system (|25|) . 

Finally, we have a dissipation term which is related to the mean temperature profile, V 2 (T). The term is nonzero 
only for the purely thermal and real modes ^"qq. In the present case this corresponds to a maximum energy of 
A$ = 0.1989%. The term is given by 

oo oo oo /*l/2 

e {T) (p x , Py , P ) = ^==y: E E a EW wSU*22U)' *o,o *fa <*•>> 



only if = p y = 0, otherwise e/n (j>x,Py,p) = 0. In Eq. (|35[) . we used again (•) = d(-)/dz. This completes the 
discussion of the different terms that arise due to the Galerkin projection of the Boussinesq equations onto the POD 
modes. 
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FIG. 2: (Color online) Reconstruction of the vertical profiles of mean convective heat flux (top) and temperature deviation 
from the linear profile (bottom). Data from DNS are compared with the two LDMs, Ml & M2, as well as with the complete 
set of POD modes (M = 15708) obtained from the snapshot analysis. 




FIG. 3: Estimation of the maximum eddy viscosity and diffusivity following the procedure of Ref. 8. (Top) Ratio Di/ei for the 
LDM model M2 with 310 POD modes. The dashed line marks the maximum of the ratio. (Bottom) Di and e; are separately 
shown for the same data. 
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Mo del 


M 


Tlx and xiy 




% Energy 


Ml 


210 


\rix \ + K| < 4 


1 < (n) < 10 


76.894 


M2 


310 


|^| + |Wy| < 5 


1 < (n) < 10 


82.202 


Total 


15708 


< n, < 16 , 
-16 < n y < 16 


1 < (n) < 28 


100 



TABLE I: Parameters of the two LDMs denoted as Ml and M2. We list the range of the horizontal wavenumbers n x and n y and 
the so-called quantum number (n). In the last two lines of the table, the total number of POD modes that has been calculated 
is given. It follows from the maximum range of horizontal wave and vertical quantum numbers. The horizontal wavenumber 
in x-direction starts from zero due to symmetry of Eq. (|16|) . The total number of POD modes from the snapshot analysis is 
thus M = 17 x 33 x 28 = 15708. 



B. Truncation to a low-dimensional model 

The integrals along the z-axis contained in the coefficients of Eq. (|25|) are evaluated on the computational grid 
of the DNS. Since the integrands are discrete functions of the wave and quantum numbers, we include all degrees 
of freedom with < n x < 16, — 16 < n y < 16, and 1 < (n) < 28. This results in a maximum number of POD 
modes of M = 15708. Out of this set of POD modes, we select small subsets of the most energetic POD modes which 
corresponds to the large-scale structures of the convection dynamics. 

The choice of M — 210 for our first LDM denoted as Ml follows from the restriction to modes with \n x \ + \n y \ < 4 
and 1 < (n) < 10 as indicated in TabJH This level of truncation builds on experiences from similar studies by Aubry 
et al£, Moehlis et al. 6 , Holmes et aZi 30 i 32 and Podvin^ 3 - in wall-bounded shear flows. They reproduced successfully 
the dynamics close to the walls for a range of values of the so-called Heisenberg parameter rj (which will be discussed 
below). Modes with n x = in streamwise and — 5 < n y < 5 in spanwise direction were taken. As explained by Holmes 
et al££, this choice is due to the fact that for the given spanwise domain length, the cross-stream interactions that 
contribute to the observed bursts of the velocity are well reproduced with at least five nonzero modes. Moreover, 
when we choose small quantum numbers, (n) < 5, the present LDM for convection yields solutions which decay 

monotonically with time. Our model and the resulting degeneracy restrictions for the average field-modes Ag did 
not allow us to take (n) < 8. A significant improvement of the LDM dynamics is obtained for (n) < 10. This is 
due to the fact that the two additional average field-modes with A 9 q = Ag 1 "" 1 have a degeneracy of 2 and will be 
incorporated. A second LDM called M2 was introduced with range of horizontal wave numbers \n x \ + \n y \ < 5 (see 
Tab.|T|. The latter LDM will be used for most of the following studies. 

Figure [2] represents the Reynolds shear stress (w9(z)) and the average temperature profile (T(z)) as reconstructed 
from the two LDMs. The calculation of both profiles is done by integration of Eqns. and ([501) . respectively and 
explained in section UlI Al As can be seen in the figure, the convergence to the vertical DNS profile for (w6(z)) is very 
slow. Only the significant enhancement of the degrees of freedom up to M = 15708 results in an excellent agreement 
with the DNS profiles, for the present Rayleigh and Prandtl numbers. 

The modes of the LDMs are however not the modes which contribute dominantly to thermal and kinetic energy 
dissipation. The missing couplings to the small-scale dissipating modes causes numerical stability problems of the 
LDM. Various methods have been proposed therefore to stabilize the truncated low-dimensional dynamical system 
as we have discussed in the introduction. Our studies showed that the model introduced by Cazemier et al& worked 
best. This model introduces a closure based on the mean energy balance as derived from Eq. (f2"5]) which can be 
rewritten as 

m q 

k p P lpy = ^2 A (Px,Py,n,p)a ( p n 2 Py + N(p x ,p y ,p) + e {T) {p x ,p y ,p) . (36) 

71=1 

where Mq is the highest quantum number in the LDM (Mq < 28) . The three linear terms in Eq. (|25[) are summarized 
to 

MPx,Py,n,p) = p(p x ,p y ,n,p) + e(p x ,p y ,n,p) + 5(p x ,p y ,n,p) . (37) 

From Eq. (|36[) . one can derive an equation for the total energy by multiplication with a p n 2 Py and summation over all 
quantum and wavenumbers. An additional linear damping term, D(p x ,p y ,p), is then quantitatively determined from 
the requirement that the mean total energy of the extended dynamical system is in a statistically stationary state, 
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y ] ' P \ m xi m yiPxiPy){ a [ n } tm y a px-mx,Py-m y a ^p],p y ) 



n .m :ni,..:}n ,. 



+ (MPx,Py,P,p) +D(Px,Py,p)\ ( a p x ,p y a p*,p y ) + £(T){Px,Py,p)(a p P J* 



(38) 



Note, that in this real equation, the last term on the right hand side is nonzero only for the purely thermal modes 
with p x —Py — 0. Also, due to orthogonality, the last two indices of A have to be equal. 

In the so-called Heisenberg dissipation model by Aubry et alA the action of the neglected modes on the ones 
contained in the LDM is represented in an average sense, namely as a function of the dynamics of these coherent 
structures. An even simpler approach was chosen by Omurtag and Sirovich.— They simply introduced a constant 
empirical viscosity coefficient for turbulent channel flows, which has a similar effect as the Heisenberg eddy viscosity 
of Aubry and co-workers. First, we will apply here the constant eddy viscosity-diffusivity, but estimate the magnitude 
of r\ by the method of Cazemier et al^ as described above. Since Pr = 0.7 and thus v sa k, we also use the same 
amplitudes of r/ for all fields. Equations (|2~5j) follow to 



Mq Mq 

' l vlvy = P(P^Py> n 'P) a p n lp y + ( 1 + r l)J2^{Px,Py,n,p)a ( p ll 



I'll 



n=l 



Mq 



N{p x ,Py,P) + ^2^(p x ,Py,n,p)a { p l\ py +€(T)(Px,Py,p) 



n=l 



after truncation and addition of an eddy viscosity and diffusivity term 

r\ = arj with afK. 
The constant eddy viscosity-diffusivity 77 is determined by 



max 



D(p x ,P y ,p) 



t(Px,p y ,n,p)5 np 



(39) 



(40) 



(41) 



Figure [3] shows the damping term Di for the POD modes of M2 and its ratio with the corresponding dissipation ti. 
In these plots i — {p x ,P y ,p) represents the POD mode i of the LDM. The modes are ordered by decreasing energy 
content (see also the second column of table [TTJ) . Note also that in Eq. (|38|) the dissipation term appears only when 
n — p. In the top panel of the figure, the modes are grouped in accordance with their degeneracy as given in Tab. [TTl 
the first two modes are (0, 1, 1) and (1, 0, 1) (the other two are obtained from complex conjugation), the second pair 
(0,1,2) and (1,0,2), etc. 

The order of magnitude of the eddy viscosity and diffusivity term is determined by the maximum value of the ratio 
Di/ei, rj « 2.698 for i — 8. It is shown in the next section that all the regimes of interest, and their typical solutions, 
can be obtained by a variation of the real prefactor a in the range a sa 0.87 — 1.0 for Ml and a « 0.87 — 1.04 for M2. 
However, as pointed out by Kalb and Deane^, this damping term Di can change sign in contrast to ti and can thus 
add as an additional production term. 

It is also evident from Fig. [3] that the Heisenberg eddy viscosity-diffusivity introduces an overwhelming damping 
for the less energetic modes. This will cause stationary convection solutions for both LDMs. These limitations of the 
model with constant eddy viscosity-diffusivity make it necessary to extend the LDM to a mode-dependent or modal 
eddy viscosity-diffusivity coefficient, -q(i). Such a model matches the decreasing ratio Di/ei. 



C. Structure of the POD modes 



Before we discuss the time evolution of the LDM with constant and modal eddy viscosity-diffusivity, we want to 
describe the structure of the POD modes. Table HIl presents the fifteen most energetic modes ^m;n x ,ri H which arc 
represented by the triplet (n x , n y ,n). Their decreasing energy content is shown in the fourth and fifth column where 

the individual eigenvalues An™^,n H and their ratio with respect to Aq 1 ^ are listed, respectively. In addition, we show 
their degeneracy and the resulting share in the total energy content. 

The POD analysis has been done for the range of horizontal wavenumbers of < n x < 16, — 16 < n y < 16. The 
order of magnitude of the resulting eigenvalue spectrum is in qualitative agreement with those from Smith et al. — ^ 
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TABLE II: The first fifteen most energetic POD modes as obtained from the snapshot analysis. The corresponding index in 
the LDM, the eigenvalues, their ratio with respect to the first eigenvalue, the degeneracy of the modes, and their percentage of 
the total mean energy content (kinetic energy plus scalar variance) are also given. Ordering is with respect to the eigenvalue 
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FIG. 4: (Color online) The wall-normal dependence of the components of the four most energetic POD modes which is given 
by (f>m}n x ,n y (z). The index m = 1,2,3,4. The particular mode is indicated in each panel. Linestyles for different components 
are the same in all four subfigures. 



and Moehlis et al. — for plane Couette flow, as well as those from Deane and Sirovich3 ? for Rayleigh-Benard convection 
in a Cartesian cell. The first few modes carry the major share in the total energy and are followed by a tail of slowly 
decaying and energetically much less significant degrees of freedom. For example, in the plane Couette flow case^ 
the first mode carries about 68% of the total energy, whereas in the RB convection case22 about 39% at Ra ps 46000. 
As shown in Refi^, the number of energetically significant modes in RB convection increases rapidly with increasing 
Rayleigh number. Thus, at Ra = 1.03 x 10 5 , our primary POD mode only represents about 31% of the total mean 
energy. The slow decay of the spectrum results in the POD mode with the triplet (0, 5, 1) which is No. i = 101 in M2 
(index i corresponds with the second column of Tab. [TTJ) still carrying 0.28% of the total energy. 
Figure [4] displays the wall-normal dependence of the components of the eigenvectors of the four most energetic 
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POD modes, which are 0^-0,1' Con 0m;i,i> an( ^ an 01 them representing a vertical circulation whose axes, 

due to their four-fold degeneracy, are along x or y, or in the case of 4>^m-i i along the horizontal diagonals of the 
periodic box, respectively. We found in our analysis that the convergence of the snapshot method is very slow. This 
causes for example the slight differences between x and <p^\ x as seen in Fig. [U Figure [5] presents the wall- normal 

functionality of further POD modes. For example, the POD mode belonging to (Pm-o o ' s a purely mechanical mode 
with degeneracy 2 and vanishing components m = 3,4. The two remaining horizontal components (to = 1,2), which 
are denoted as 0^-0 o an d <^m-o 0' sa tisfy orthogonality (see panel (a)). 

Another interesting mode is the sixth most energetic POD mode belonging to (/> m 'o i- ^ represents a pumping 
motion in the x or y directions and is the hrst of the most energetic POD modes with a non-vanishing vertical 
vorticity component^ The mode is shown in panel (b) of Fig. [5] Furthermore, we also display a purely thermal 
mode ^m-o o with degeneracy 1 in panel (c) and the mode 4&. 8 which shows the growing influence of the viscous 
and thermal boundary layers close to the walls at higher horizontal wave numbers in panel (d) of the figure. 
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z z 

FIG. 5: (Color online) Wall-normal dependence of further POD modes. The legend in panels (c) and (d) is identical to panel 
(b). 

Figure |5] displays the spatial structure of the velocity field corresponding to the functions 1 and 1 ■ They 
are obtained by summing all individual eigenfunctions as given by Eq. (|12p over their degeneracy, i.e., modes that 
have the same energy content. Figure [7] illustrates the isosurfaces of the temperature fluctuations 0(x, t) that are 
captured by the fourth component of the same POD modes. We show 1 in the top panel and $ 4 X ] 1 in the bottom 
panel of the figure. 

One important aspect of the POD mode analysis and LDM setup is the question of how well the largest structures 
and their dynamics can be represented. Here, we display the reconstruction of a DNS snapshot by the mode set. In 
Fig. [8j we compare the velocity field as reconstructed with the POD modes $m;n»,n v from model M2. The top panel 
in the figure shows the reconstruction and the bottom panel the full DNS snapshot. The same analysis is repeated for 
the temperature fluctuations at the same instant. Figure [3] shows the temperature field again reconstructed with M2. 
Although, not all details are reproduced, both figures indicate that the most important structures of temperature 
and velocity are captured by M2. Our analysis showed that the temperature fluctuations converged slower than the 
velocity field, ft is known that the temperature field forms so called thermal plumes - fragments of the thermal 
boundary layer that detach from the cooling and heating plates and move into the bulk (see e.g. Refs. 35, 36 and 
37) . These fine-scale filamented structures carry the heat across the cell. 
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FIG. 6: (Color online) Vector plots of the first and third most energetic POD modes 5^. 1 (top) and 1 (bottom). The 
velocity field is therefore projected into three y-z planes. 




X 



FIG. 7: (Color online) Isosurfaces of the first and third most energetic POD modes $^ ^ (top) and ^ (bottom) which 
describe the temperature fluctuations. The isosurfaces are taken at the level 0.01 9 C . In addition, we plot contours at the 
sideplanes. Red color stands for 1.7 6 C and blue for -1.7 6 C . 
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FIG. 8: (Color online) A snapshot of the velocity field u(x, to) obtained from the DNS (bottom) is reconstructed with the POD 
modes of low-dimensional model M2 (top). 



D. Time evolution of LDM with constant eddy viscosity— diffusivity 

In the following, we will discuss the evolution of the LDM. This section is for the case with the constant eddy 
viscosity-diffusivity, i.e. 77 = arj where rj is given by (1411) . and a 6 K. 

As was found in Ref. 8 for the long-term integration of their LDM, this closure can lead to a statistical equilibrium 
state which accumulates too much energy. One can overcome this behavior in parts following a work by Kalb and 
Deane^, who kept the original gradient of the mean temperature field, , in the evaluation of the last two terms of 
Eq. ([25jl . Consequently, the equations for G and em become 

%{ Px ,p y ,n,p) = -u c J ^ ^-^L^^p^Pv dz - ( 42 ) 

and, 

/1/2 J2(T) 
-Tj-C dz ■ ( 43 ) 
-1/2 az 

The long-time evolution of the LDM requires the time integration of the ODE system A fourth-order Runge- 

Kutta scheme is applied. Our studies found that for a ig 0.87 the ODE systems for Ml and M2 become unstable 
which is triggered by the most energetic modes that accumulate energy which cannot be transferred sufficiently fast to 
small scales. For a > 0.87, we still obtain a regime of the LDM carrying too much energy. With increasing prefactor 
a, this energy surplus at the first POD modes decreases up to a threshold which enforces the whole dynamical system 
into a stationary state. For Ml this sets in at a > 1.0 and for M2 at a > 1.04. In Fig. [TU1 we compare therefore the 
total energy of both LDMs for different values of a with the original DNS data. After a relaxation phase both models 
reach a statistically stationary regime. While model Ml always yields energy time series below that of the original 
DNS, model M2 comes closer to the evolution of the total energy from the original DNS. Note that the fluctuations 
of the total energy in both cases are significantly larger than for the original DNS data. This is due to the fact that 
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FIG. 9: (Color online) Temperature fluctuation field 8(x,t) of a DNS snapshot (bottom panel) is reconstructed with the POD 
modes of low-dimensional model M2 (top panel). The isosurfaces are taken at the level 0.4 9 C . In addition, we plot contours at 
the sideplanes. Red color stands for 1.7 9 C and blue for -1.7 6 C . 

only the largest-scale modes are kept in the model and coupled with each other. The exchange of energy among them 
can cause larger variations. The figure clearly indicates that Ml falls short in representing the long-term dynamics of 
the convective flow. In model M2 however it is possible to obtain a total energy in the range of the original data. It 
can be concluded that a further reduction of degrees of freedom below that of Ml is thus not possible. 

Figure Qj] shows the long-time behavior of \y affi^y an d &^2) f° r ^-2 in the stationary overdamped state. 

This fixed point of the LDM is connected with a stationary pattern of the velocity and temperature fields. It is now 
clear how a constant eddy viscosity-diffusivity model produces the observed behavior. The second mode (dashed line) 

is accumulating too much energy )l l a (o\)l)' smce it is n °t damped strongly enough. Similar behavior of the 
LDM has been found and analyzed in detail by Aubry et al£>. They detected a similar fixed point for a specific range 
of their constant eddy viscosity. 

Since the overdamped stationary state is not of interest for our study, the values of a have to be chosen smaller than 
these limits. Figure [T2l shows the time evolution of the modal amplitude for the most energetic POD mode, %\(t), 
obtained from the LDMs Ml and M2, respectively at a smaller value of a. Data are again compared with DNS time 
series which are obtained by projection of the snapshots on the particular modes. Clearly, at a = 0.93 the agreement 
of the real part is very good for M2, even for long term evolution. The amplitude of the imaginary part in M2 is 
too large, the one in Ml is comparable with the DNS. While the amplitudes partly agree, the temporal behavior of 
the modes differs qualitatively. As observed from Fig. [T^J the real and imaginary parts of the expansion coefficients 
vary in a limited range (in parts periodically) once the initial relaxation to a statistically stationary state is finished. 
This is in contrast to the projection of the DNS snapshots on the modes. It indicates that the additional dissipation 
has a strong impact on the dynamics of the large-scale degrees of freedom. The figure unravels the shortcoming of 
the present straightforward and simplest closure: all modes are affected by the same additional dissipation, the ones 
that have many couplings within the model as well as those with much less mode interactions. The constant eddy 
viscosity-diffusivity establishes an additional flux from large resolved to small unresolved scales and seems not to 
allow a back-scatter which is important and known from other (subgrid-scale) closures^ 
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FIG. 10: Time series of the total energy for LDM Ml (top panel) and M2 (bottom panel). Different magnitudes of the parameter 
a in the constant Heisenberg eddy viscosity-diffusivity are applied in both models and given in the legend of the bottom panel. 
Note that a = 1.04 is not taken for Ml. 
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FIG. 12: Time series for the real (top) and imaginary (bottom) parts of of the two LDMs. The parameter a is set to a 
value of 0.93. 



E. Time evolution of the LDM with modal eddy viscosity— diffusivity 




FIG. 13: Modal eddy viscosity-diffusivity following the maximum values of Di/ti. The real prefactor is now slightly larger 
and given by a — 1.146 for the solid line. For comparison, we add the constant value of rj from Fig. [3] Model M3 is used for 
evaluation of M2 only. The inset shows the degeneracy restrictions that have to be included for the first 41 modes. 



A refinement of the constant eddy viscosity-diffusivity model, that overcomes the shortcomings from above, is 
possible when switching to the modal eddy viscosity-diffusivity. Figure [13] shows the maximum values of the ratio 
Di/a that were already presented in the top panel of Fig. [3] anew for the eddy viscosity-diffusivity of model M2. 
The sudden decrease in the maxima for i w 140 is a consequence of the truncation. We confirmed this after plotting 
the modes resulting from a slightly larger model -denoted as M3- that contains 430 modes with \n x \ + \n y \ < 6, 
1 < (n) < 10 (crossed symbols in the figure). The solid line is a fit which is given by 

= a/3( 7 r • (44) 
We kept the variable real prefactor a in order to compare this case with the former constant eddy viscosity-diffusivity. 
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Clearly the new function /3(7) z is accounting now for the constant rj in the Heisenberg dissipation model. Data fit 
then well for /3 = 2.639, and 7 = 0.99372, in the range for i < 140. It turned out to be necessary to add some more 
damping through the factor a > 1 and to check the results for the ensemble average and the transient total energy. 
In agreement with Ref. 8, it yields a positive definite dissipative damping term Di for all modes. The inset shows the 
degeneracy for the first 41 modes. 
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FIG. 14: Top panel: Mode spectra obtained from a run of M2 with a = 1.146 and the original POD of the DNS data. 
Xi — (a^ Vy a? * Py ) with i — (px,p y ,p) as in section UlIBI Bottom panel: Time series of the total energies. 

Figure [14] (top panel) displays the energy spectrum of model M2 after integration of the ODE system following 
the modal dependent closure model (solid line). It is compared with the original POD spectrum of the DNS-data 
(dashed line). Due to the high accuracy at the tail, the first 120 modes are shown only. It is obvious that a further 
improvement in the accuracy of the M2-spectrum can be attained considering the effect of the C/-symmetry group (see 
Appendix [A} on the time evolution coefficients, as it was done for the full DNS data. The analytical determination 
of the effect of these symmetries, as was studied by Smith et al&l for turbulent plane Couette flow, is beyond the 
scope of the present study and must be addressed in the future. Differences appear for the two most energetic modes 
(0,1,1), (1,0,1) which have slightly smaller energy. Overshoots for modes No. 12 with (n x ,n y ,n) — (0,1,3) , 25 with 
(1,-1,3), 26 with (1,1,3), and 29 with (0,0,3) are found. Undershoots are detected for modes No. 9 with (0,0,1), 10 
with (0,0,2), 17 with (1,0,4) through 20 with (2,0,2), and 30 with (1,1,4). 

The bottom panel, which is analogous to the bottom panel of Fig. [101 shows the instantaneous total energy for 
500 time units. The agreement between DNS and M2 is now significantly better. Energy remains largely fluctuating 
and yields an ensemble average equivalent to 91.11% of the amplitude of the DNS, 5.58% smaller than the 96.19% 
corresponding to constant eddy viscosity-diffusivity at a — 1.04 (bottom panel of Fig. [TO)) . 

In Fig. Qjj] we show the turbulence statistics as obtained from a long-time run of model M2 with modal eddy 
viscosity-diffusivity. Plane and time averaged vertical profiles of (it 2 ), (v 2 ), (w 2 ), and (6* 2 ), computed following 
equations analog to (|33|) . Time average is taken over the first 500 time units. In each case, all the main features of the 
original DNS-profiles are reproduced, and the truncation yields a reasonable accuracy. The profiles are reproduced 
qualitatively well. 

Figure [TBI shows the time evolution of the modal amplitude for the most energetic POD mode, a^\(t), obtained from 
M2 using the modal eddy viscosity-diffusivity. Comparison with the DNS time series shows a reasonable agreement 
of the real and the imaginary parts until about t = 90. Furthermore, the whole temporal evolution of the coefficients 
is now much closer to those of the DNS. In contrast with Fig. I12[ now the real and imaginary components vary in 
ranges as wide as the extent of variation for the projection of the DNS snapshots on the modes. The relaxation into 
a quasi-periodic time variation as in the case with constant eddy viscosity-diffusivity as now absent. 
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FIG. 15: Turbulence statistics (u 2 ), (v 2 ), (w 2 ), and {O 2 } versus z as computed from the time evolution of model M2 and 
compared with the DNS. Linestyles are the same in the four subfigures. The profiles are obtained by taking plane-time 
averages. 




Figure [TTl shows the total temperature field, 0(x, t) (see Eq. (HJ)), reconstructed with M2, at the two instants of 
time t = 100 and 400. We clearly identify the mushroom-shaped isosurfaces, which belong to the same type as the 
ones presented in Fig. [9] for the fluctuations, and were previously observed (see e.g. Refs . 36 i ) for RB convection in 
cylindrical containers for Ra = 10 5 -10 9 . Velocity snapshots at the four instants t = 100, 200, 300 and 400 are shown 
in Fig. 1181 It is observed that the flow is characterized by strong up- and downward flows, which are in line with the 
enhanced fluctuations of the total energy in Fig. [14] 
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FIG. 17: (Color online) Isosurfaces of the total temperature field O at the two instants of time: t = 100 (top row) and t = 400 
(bottom row). Isosurfaces at O(x.,t)/0 C = 0.7 (left column with top view) and —0.7 (right column with a view from below). 
Data are obtained from model M2. 




FIG. 18: (Color online) Velocity field snapshots at four instants of time: (a) t = 100, (b) t = 200, (c) t = 300, (d) t = 400. 
Data are obtained from model M2. 
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IV. SUMMARY AND OUTLOOK 



We have studied a low-dimensional model of turbulent Rayleigh-Benard convection in a Cartesian slab. The POD 
modes which form the basis of our model have been obtained by a so-called snapshot method from a record of 320 
statistically independent realizations of a DNS of convective turbulence for the same geometry. Temperature and 
velocity field fluctuations have to be considered therefore as a common four-component vector field where the mode 
selection is done with respect to the total energy in the convective flow, i.e. kinetic energy plus thermal variance. The 
Navier-Stokes-Boussinesq equations are then projected onto the POD modes. The Galerkin projection is truncated 
at two different levels and results in the low-dimensional models denoted as Ml and M2. 

Our results can be summarized as follows. The LDMs have to be stabilized by an additional eddy viscosity-diffusivity 
r\ that assures that the generated energy can be dissipated since the small-scale degrees of freedom are missing in the 
model. This observation is in line with existing works on flows in channels or cavities. First, we introduced a constant 
r\ whose order of magnitude, 77, was determined as the global maximum quotient of the damping term, Di, and the 
diffusivity The calculation of term Di is based on the requirement to have a statistically stationary dynamics in 
the LDM. We have then studied the dynamics of our LDM as a function of the level of truncation and the magnitude 
of rj. Similar to the works in simple wall-bounded flow or cavities, we observe a convergence into a stationary regime 
for amplitudes of rj that are too large. This regime can be considered as a fixed point which is however not in the 
focus of the present study. We showed also that the effect of r\ on the dynamics of the largest-scale modes is to force 
them into a state with small fluctuations after the passage of a longer transient. 

Alternatively, we introduced a modal eddy viscosity-diffusivity, rj(i), by fitting the algebraic power law of the form 
P('y) 1 to all the local maxima of the quotient Di/ei, and considering the restrictions due to the degeneracy of the most 
energetic modes. For the model M2, the long time integration of the ODE system yields solutions with remarkable 
accuracy in the value of the ensemble average of the total energy (5.6% below the value for constant rj). Also, the 
energy spectrum of the POD, is very well reproduced, especially in the tail for modes with higher indexes. The 
vertical profiles of plane-time averaged fluctuations agree qualitatively with those from DNS. Characteristic coherent 
structures of convection, such as thermal plumes, are reproduced. We can thus conclude that the second approach 
with the modal eddy viscosity-diffusivity can model the long-term dynamics of turbulent convection qualitatively 
well. We wish to stress here again that this was in the focus of the present work, namely how far we can advance 
with a least set of POD modes. 

The question is interesting and important in view to more complex situations, e.g. the problem of mixed convection 
in indoor ventilation systems. Can the same POD framework (with a mathematical foundation) be carried over to 
more complex convection flows? A big advantage of the present turbulent convection case in the Cartesian box with 
periodic side walls is that we have 16 symmetries that significantly enhance the data base. In view of applications 
in more complex geometries, this indicates that a similar approach might be much more complicated and could 
enhance the limitations that showed up already for the present turbulent flow. Nevertheless, since the questions, for 
example with the long-term behavior of the large-scale circulation in turbulent convection^ or indoor ventilation, are 
important, we believe that it is still interesting to further follow this route of LDM development based on the POD 
framework. Some of these efforts will be hopefully presented in the near future. 
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The numerical simulation of Rayleigh-Benard convection for the case of a square section, L x — L y , generates a 
maximal amount of symmetry. These discrete symmetries form a group G of eight elements^ 
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Appendix A: Symmetry considerations 



G = {/, R, R 2 , R 3 , F, FR, FR 2 , FR 3 } 



(Al) 



whose generators 



are the rotation by 90° 



R(x, y, z, u, v, w, 9) = (— y, x, z, —v, u, w, 9) 



(A2) 
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and the reflection in x 

F(x, y, z, u, v, w, 6) = (-x, y, z, -u, v, w, 6) (A3) 
Furthermore, another symmetry group {I, Z} acts on the vertical direction, where Z is the reflection in z 

Z(x, y, z, u, v, w, 9) = (x, y, -z, u, v, -w, -6) (A4) 
resulting in a symmetry group of sixteen elements 

£ = {G,ZG} (A5) 

when combined. Since each element of the symmetry group generates a possible flow, the ensemble is enlarged by a 
factor of sixteen, thus increasing the accuracy of any statistical evaluation of the flow. 
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